Estimating Seasonal Behavior States from Biologging Sensor Data

data(kotzeb0912_locs)

locs <- dplyr::filter(kotzeb0912_locs,instr=="Mk10") %>% 
  dplyr::arrange(deployid,unique_posix)
# speedfilter using the paralell package for multi-core speed
cfilter<-mclapply(split(locs,locs$deployid),function(x) sdafilter(
  lat=x$latitude, lon=x$longitude, dtime=x$unique_posix,
  lc=x$quality, ang=-1,vmax=5),mc.preschedule=F,mc.cores=3)
cfilter<-do.call("c",cfilter)
cfilter<-as.vector(cfilter)
locs$filtered <- cfilter

data <- filter(locs,filtered=="not", !is.na(error_semimajor_axis)) %>% arrange(deployid,unique_posix) %>% as.data.frame(.)

data(kotzeb0912_gps)

kotzeb0912_gps <- kotzeb0912_gps %>% 
  select(deployid,date_time,latitude,longitude) %>% 
  filter(!is.na(latitude)) %>% 
  rename(unique_posix=date_time) %>% 
  mutate(error_semimajor_axis = 50,error_semiminor_axis = 50, error_ellipse_orientation = 1) 

data <- data %>% bind_rows(kotzeb0912_gps) %>% 
  arrange(deployid,unique_posix)

data <- as.data.frame(data)
coordinates(data) = ~longitude+latitude
proj4string(data) = CRS("+proj=longlat +datum=WGS84")

data=spTransform(data, CRS("+init=epsg:3571"))

## Loop over PTT and fit CTCRW models
ids = unique(data@data$deployid)

library(doParallel)
n.cores <- detectCores()
registerDoParallel(cores=3)

model_fits <-
  foreach(i = 1:length(ids)) %dopar% {
  id_data = subset(data,deployid==ids[i])
  diag_data = model.matrix(
      ~ error_semimajor_axis + error_semiminor_axis + error_ellipse_orientation,
      id_data@data
    )[,-1]
    
    id_data@data = cbind(id_data@data, 
                         argosDiag2Cov(
                           diag_data[,1], 
                           diag_data[,2], 
                           diag_data[,3]))
    
    init = list(a = c(coordinates(id_data)[1,1],0,
                      coordinates(id_data)[1,2],0),
                P = diag(c(5000 ^ 2,10 * 3600 ^ 2, 
                           5000 ^ 2, 10 * 3600 ^ 2)))
    
  fit <- crwMLE(
      mov.model =  ~ 1,
      err.model = list(
        x =  ~ ln.sd.x - 1, 
        y =  ~ ln.sd.y - 1, 
        rho =  ~ error.corr
      ),
      data = id_data,
      Time.name = "unique_posix",
      initial.state = init,
      fixPar = c(1,1,NA,NA),
      theta = c(log(10), 3),
      initialSANN = list(maxit = 2500),
      control = list(REPORT = 10, trace = 1)
    )
    fit
  }

names(model_fits) <- ids

print(model_fits)

predData <- foreach(i = 1:length(model_fits)) %dopar% {

  model_fits[[i]]$data$unique_posix <- lubridate::with_tz(
    model_fits[[i]]$data$unique_posix,"GMT")
  predTimes <- seq(
    lubridate::ceiling_date(min(model_fits[[i]]$data$unique_posix),"day"),
    lubridate::floor_date(max(model_fits[[i]]$data$unique_posix),"day"),
    "6 hours")
  tmp = crwPredict(model_fits[[i]], predTime=predTimes)
}


predData <- dplyr::bind_rows(predData) %>% 
  filter(locType=="p")

predData$predTimes <- intToPOSIX(predData$TimeNum)

data(kotzeb0912_timelines)

kotzeb0912_timelines <- kotzeb0912_timelines %>% 
  filter(deployid %in% unique(predData$deployid))

ids <- unique(kotzeb0912_timelines$deployid)
timelineData <- foreach(i = 1:length(ids)) %dopar% {
  t_dat <- kotzeb0912_timelines %>% 
    filter(deployid == ids[i]) %>% 
    arrange(deployid,datadatetime)
  t_seq <- data.frame(datadatetime = seq(
    lubridate::ceiling_date(min(t_dat$datadatetime),"day"),
    lubridate::floor_date(max(t_dat$datadatetime),"day"),
    "1 hour"),deployid=ids[i])
  t_dat <- merge(t_dat,t_seq,all=TRUE)
  rep_len <- ceiling((nrow(t_dat)-3)/6)
  t_group <- c(rep(1,3),rep(2:rep_len,each=6))
  t_group <- data.frame(group=t_group[1:nrow(t_dat)])
  t_dat <-cbind(t_dat,t_group) %>% filter(!is.na(group))
}

timelineData <- dplyr::bind_rows(timelineData) %>% 
  group_by(deployid,group) %>% 
  summarise(percent_dry=mean(percent_dry,na.rm=TRUE),
            datadatetime=ceiling_date(median(datadatetime),"hour")) %>% 
  select(deployid,datadatetime,percent_dry)

predData <- dplyr::inner_join(predData,
                              timelineData,
                              by = c("predTimes" = "datadatetime", 
                                     "deployid")) %>%
  tbl_df() %>%
  select(deployid,predTimes,mu.x,mu.y,percent_dry) %>% 
  arrange(deployid, predTimes) %>% 
  rename(datadatetime=predTimes)

data(kotzeb0912_depths)

diveData <- kotzeb0912_depths %>% 
  filter(limits != "10.000000") %>% 
  filter(num_dives < 80) %>% 
  mutate(datadatetime = datadatetime + lubridate::hours(3)) %>% 
  group_by(deployid,datadatetime) %>% 
  summarise(num_dives=sum(num_dives,na.rm=TRUE))

predData <- dplyr::left_join(predData,
                             diveData) %>% 
  tbl_df() %>% 
  arrange(deployid,datadatetime)

predData <- predData %>% 
  group_by(deployid) %>% 
  mutate(x_disp = order_by(datadatetime, mu.x - lag(mu.x)), 
         y_disp = order_by(datadatetime, mu.y - lag(mu.y))
  ) %>% 
  filter(!is.na(x_disp) | !is.na(y_disp))

anglefun <- function(xx,yy,bearing=TRUE,as.deg=FALSE){
  ## calculates the compass bearing of the line between two points
  ## xx and yy are the differences in x and y coordinates between two points
  ## Options:
  ## bearing = FALSE returns +/- pi instead of 0:2*pi
  ## as.deg = TRUE returns degrees instead of radians
  c = 1
  if (as.deg){
    c = 180/pi
  }
  
  b<-sign(xx)
  b[b==0]<-1  #corrects for the fact that sign(0) == 0
  tempangle = b*(yy<0)*pi+atan(xx/yy)
  if(bearing){
    #return a compass bearing 0 to 2pi
    #if bearing==FALSE then a heading (+/- pi) is returned
    tempangle[tempangle<0]<-tempangle[tempangle<0]+2*pi
  }
  return(tempangle*c)
}

predData <- predData %>% 
  mutate(bearing = anglefun(x_disp,y_disp,as.deg=TRUE))

predData %>% group_by(deployid) %>% 
  filter(row_number(datadatetime) == 1)

predData %>% group_by(deployid) %>% 
  filter(row_number(datadatetime) == max(row_number(datadatetime)))

predData <- predData %>% 
  mutate(deploy_day = ifelse(lubridate::yday(datadatetime)<150, 
                      lubridate::yday(datadatetime)+365 +
                        lubridate::hour(datadatetime)/24,
                      lubridate::yday(datadatetime) +
                        lubridate::hour(datadatetime)/24)) %>% 
  mutate(corrupt = ifelse(deployid == "EB2009_3001_06A1332" & 
                            datadatetime > lubridate::ymd_hms("2010-01-28 03:00:00"),
                          TRUE,FALSE
                          )) %>% 
  filter(abs(y_disp) < 100000, abs(x_disp) < 100000) %>% 
  filter(!corrupt)

save(predData,file='predData.rda')

library(ggplot2)
library(uswebr)

p1 <- ggplot(predData, aes(x=deploy_day, y=percent_dry)) + geom_point(alpha=0.2) + facet_grid(deployid ~ .) + 
  usweb_theme()
p1

p2 <- ggplot(predData, aes(x=deploy_day, y=num_dives)) + 
  geom_point(alpha=0.2) + facet_grid(deployid ~ .) + 
  usweb_theme()
p2

predData <- predData %>% mutate(direction = ifelse(
  y_disp >= 0, "north","south"))

ggplot(predData,aes(x=deploy_day, y=y_disp/1000)) +
  geom_bar(aes(color=factor(direction)),alpha=0.2,stat="identity") +
  facet_grid(deployid ~ .) + usweb_theme()